/**
 *  @file
 *  @brief   Tetra4单元
 *  @author  kangrongyao
 *  @date    2024.04.18
 *  @version v1.0
 *  @par Copyright(c):  ShenGong  Simulation corporation
 */
#include "Tetra4ElementCalculator.h"
#include <algorithm>
#include "Utility/Logging/SGLogger.h"
#include "Utility/Algebra/Vector3D.h"
#include "Utility/Algebra/MathTools.h"
#include "Utility/Algebra/Tools.h"
#include "DataStructure/Common/Id.h"
#include "DataStructure/Input/Solver/DOFDirections.h"
#include "DataStructure/Input/Solver/TenserTranform.h"
#include "DBManager/IService/Input/Node/INodeService.h"
#include "DBManager/IService/Input/Element/IElementService.h"
#include "DBManager/IService/Input/Material/IMaterialSummaryService.h"
#include "DBManager/IService/Input/Property/IPSOLIDService.h"
#include "DBManager/IService/Result/ElementalStrain/ITETRAStrainService.h"
#include "DBManager/IService/Result/ElementalStress/ITETRAStressService.h"
#include "SGFem/Pre/Assembly/Assembly.h"
#include "SGFem/Pre/Assembly/CoordinateProcessor.h"
#include "SGFem/Pre/Assembly/ModelInfoQuery.h"
#include "SolidMaterialParser.h"

using SG::Algebra::Matrixd;
using SG::Algebra::Vector3D;
using SG::Algebra::PointArray;
using SG::DataStructure::Common::Id;
using SG::DataStructure::Common::IdArray;
using SG::DataStructure::Common::Real;
using namespace SG::DataStructure::FEM;
using SG::DataStructure::Post::TETRAStrainData;
using SG::DataStructure::Post::TETRAStressData;
using namespace SG::DBManager;
using SG::FEMSolver::AssemblyEigen;
using SG::FEMSolver::CoordinateProcessor;
using SG::FEMSolver::ModelInfoQuery;
using SG::FEMSolver::Element::SolidMaterialParser;
using SG::FEMSolver::Element::Tetra4ElementCalculator;

void Tetra4ElementCalculator::Initialize (SG::DataStructure::FEM::ElementTypes type)
{
    m_type = type;
    m_shpFunc.Compute (m_integCalc.GetPoints ());
}

void Tetra4ElementCalculator::CoordTrans
(
const SG::Algebra::PointArray& pointSet,
_OUT SG::Algebra::Matrixd&     T,
_OUT SG::Algebra::Matrixd&     localCoord
)
{
    // 节点坐标矩阵化
    Matrixd basicCoord (3, 4);
    for (std::size_t iNode{0}; iNode < 4; ++iNode)
    {
        basicCoord (0, iNode) = pointSet[iNode].m_x;
        basicCoord (1, iNode) = pointSet[iNode].m_y;
        basicCoord (2, iNode) = pointSet[iNode].m_z;
    }

    // 计算单元坐标系参考基矢量 R、S、T（非标准正交），以列向量形式组成矩阵
    const Matrixd rstTrans (3, 4, {-1.0, -1.0,  1.0,  1.0,
                                   -1.0,  1.0, -1.0,  1.0,
                                   -1.0,  1.0,  1.0, -1.0});
    const Matrixd rst = basicCoord * rstTrans.Transpose ();

    // 参考基矢量标准正交化，即得单元坐标系变换矩阵
    T = SG::Algebra::orthonormalize (rst);

    // 单元节点变换到以节点 1 为原点的单元坐标系下
    localCoord = Matrixd (3, 4);
    for (std::size_t iNode{1}; iNode < 4; ++iNode)
    {
        localCoord (0, iNode) = basicCoord (0, iNode) - basicCoord (0, 0);
        localCoord (1, iNode) = basicCoord (1, iNode) - basicCoord (1, 0);
        localCoord (2, iNode) = basicCoord (2, iNode) - basicCoord (2, 0);
    }
    localCoord = T.Transpose () * localCoord;
}

void Tetra4ElementCalculator::ComputeInternalCoordTrans (const SG::Algebra::PointArray& pointSet, _OUT SG::Algebra::Matrixd& TS)
{
    Vector3D xVec{(pointSet[1].m_x + pointSet[2].m_x + pointSet[3].m_x) / 3.0 - pointSet[0].m_x,
                  (pointSet[1].m_y + pointSet[2].m_y + pointSet[3].m_y) / 3.0 - pointSet[0].m_y,
                  (pointSet[1].m_z + pointSet[2].m_z + pointSet[3].m_z) / 3.0 - pointSet[0].m_z};
    Real L;
    xVec.Normalize (L);
    const Vector3D orienVector{pointSet[2].m_x - pointSet[1].m_x,
                               pointSet[2].m_y - pointSet[1].m_y,
                               pointSet[2].m_z - pointSet[1].m_z};
    SG::Algebra::computeLocalCoordTrans (xVec, orienVector, _OUT TS);
}

void Tetra4ElementCalculator::CalcBMatrix (const SG::Algebra::Matrixd& localCoord, _OUT SG::Algebra::Matrixd& B,  _OUT SG::DataStructure::Common::Real& V)
{
    static const Matrixd dNdxi (3, 4, {-1.0,  1.0,  0.0,  0.0,
                                       -1.0,  0.0,  1.0,  0.0,
                                       -1.0,  0.0,  0.0,  1.0});

    const auto dxdxi = localCoord * dNdxi.Transpose ();

    Real       detJ;
    const auto J    = SG::Algebra::inverse3X3 (dxdxi, detJ);
    const auto dNdx = dNdxi.Transpose () * J;
    V = detJ / 6.0;

    // 使用形函数的导数来构建 B 矩阵
    B = Matrixd (6, 12);
    for (std::size_t iNode{0}; iNode < 4; ++iNode)
    {
        const auto i = iNode * 3;
        B (0, i)     = dNdx (iNode, 0);
        B (1, i + 1) = dNdx (iNode, 1);
        B (2, i + 2) = dNdx (iNode, 2);
        B (3, i)     = dNdx (iNode, 1);
        B (3, i + 1) = dNdx (iNode, 0);
        B (4, i + 1) = dNdx (iNode, 2);
        B (4, i + 2) = dNdx (iNode, 1);
        B (5, i)     = dNdx (iNode, 2);
        B (5, i + 2) = dNdx (iNode, 0);
    }
}

void Tetra4ElementCalculator::ComputeStiffness
(
const SG::Algebra::Matrixd& localCoord,
const SG::Algebra::Matrixd& D,
_OUT SG::Algebra::Matrixd&  kl
)
{
    Matrixd B;
    Real V;
    CalcBMatrix (localCoord, _OUT B, _OUT V);

    kl = B.Transpose () * D * B * V;
}

void Tetra4ElementCalculator::ComputeLumpedMass
(
const SG::Algebra::Matrixd&            localCoord,
const SG::DataStructure::Common::Real& rho,
_OUT SG::Algebra::Matrixd&             lumpedM
)
{
    const auto v  = SG::Algebra::calcTetraVolumn (localCoord.Transpose ());
    const auto mi = v * rho * 0.25;
    for (std::size_t i{0}; i < 12; ++i)
    {
        lumpedM (i, i) = mi;
    }
}

void Tetra4ElementCalculator::ComputeStressStiffness (const Matrixd& localCoord, const Matrixd& eleEleCoordStress, _OUT Matrixd& kSigma)
{
    //计算单元坐标系下应力刚度矩阵
    static const Matrixd dNdxi (3, 4, {-1.0,  1.0,  0.0,  0.0,
                                       -1.0,  0.0,  1.0,  0.0,
                                       -1.0,  0.0,  0.0,  1.0});

    const auto dxdxi = localCoord * dNdxi.Transpose ();

    Matrixd G (9, 12);
    Matrixd S (9, 9);

    Real       detJ;
    const auto J    = SG::Algebra::inverse3X3 (dxdxi, detJ);
    const auto dNdx = dNdxi.Transpose () * J;
    const Real V    = detJ / 6.0;

    for (std::size_t iNode{0}; iNode < 4; ++iNode)
    {
        const auto i = iNode * 3;
        G (0, i) = dNdx (iNode, 0);
        G (1, i) = dNdx (iNode, 1);
        G (2, i) = dNdx (iNode, 2);
        G (3, i + 1) = dNdx (iNode, 0);
        G (4, i + 1) = dNdx (iNode, 1);
        G (5, i + 1) = dNdx (iNode, 2);
        G (6, i + 2) = dNdx (iNode, 0);
        G (7, i + 2) = dNdx (iNode, 1);
        G (8, i + 2) = dNdx (iNode, 2);
    }

    for (std::size_t i{0}; i < 3; ++i)
    {
        S (3 * i,     3 * i + 1) = eleEleCoordStress (3, 0);
        S (3 * i,     3 * i + 2) = eleEleCoordStress (4, 0);
        S (3 * i + 1, 3 * i + 2) = eleEleCoordStress (5, 0);
        S (3 * i + 1, 3 * i)     = S (3 * i,     3 * i + 1);
        S (3 * i + 2, 3 * i)     = S (3 * i,     3 * i + 2);
        S (3 * i + 2, 3 * i + 1) = S (3 * i + 1, 3 * i + 2);

        for (std::size_t j{0}; j < 3 ; ++j)
        {
           S (3 * i + j, 3 * i + j) = eleEleCoordStress (j, 0);
        }
    }
    kSigma = G.Transpose () * S * G * V;
}

SG::Algebra::Matrixd Tetra4ElementCalculator::ComputeCoupledMass
(
const SG::Algebra::Matrixd&            localCoord,
const SG::DataStructure::Common::Real& rho
)
{
    Matrixd coupledM (12, 12);
    Matrixd jacobianM (3, 3);
    const auto nGauss = m_integCalc.GetPointSize ();
    for (std::size_t iGauss{0}; iGauss < nGauss; ++iGauss)
    {
        const auto& dNdxi   = m_shpFunc.GetShapeFunctionDerivative (iGauss);
        jacobianM           = dNdxi * localCoord.Transpose ();
        const auto rhoDetJW = SG::Algebra::det3X3 (jacobianM) * rho * m_integCalc.GetPoint (iGauss).m_weight / 6.0;
        const auto& Ni      = m_shpFunc.GetShapeFunction (iGauss);
        for (std::size_t j{0}; j < 4; ++j)
        {
            for (std::size_t k{0}; k < 4; ++k)
            {
                const auto NiNj = Ni (j, 0) * Ni (k, 0) * rhoDetJW;
                coupledM (j * 3, k * 3)         += NiNj;
                coupledM (j * 3 + 1, k * 3 + 1) += NiNj;
                coupledM (j * 3 + 2, k * 3 + 2) += NiNj;
            }
        }
    }

    return coupledM;
}

void Tetra4ElementCalculator::ComputeElementMatries (const SG::DBManager::DBServiceFactorySP dbServiceFactory)
{
    auto&       assembly       = AssemblyEigen::GetInstance ();
    const auto& coordProcessor = CoordinateProcessor::GetInstance ();

    static const DofSet dofSet{DOFDirections::DOF_Ux, DOFDirections::DOF_Uy, DOFDirections::DOF_Uz};

    CTETRA4Data              element;
    std::vector<std::size_t> nodeIndexSet (4);
    std::vector<int>         dofGlobalIndexSet (12);
    PointArray               pointSet (4);
    IdArray                  dispCoordIdSet (4);
    std::vector<Matrixd>     TGSet (8);
    PSOLIDData               property;

    Matrixd T (3, 3);
    Matrixd localCoord (3, 4);
    Matrixd TS (3, 3);
    Matrixd TMat (3, 3);
    Real    rho;
    Matrixd ke (12, 12);
    Matrixd mass (12, 12);

    const auto& nodeService            = dbServiceFactory->get<INodeService> ();
    const auto& eleService             = dbServiceFactory->get<IElementService> ();
    const auto& propertyService        = dbServiceFactory->get<IPSOLIDService> ();
    const auto& materialSummaryService = dbServiceFactory->get<IMaterialSummaryService> ();
    const auto  eleSize                = eleService->GetElementSize (m_type);

    for (std::size_t iEle{0}; iEle < eleSize; ++iEle)
    {
        eleService->FindByIndex (iEle, _OUT element);

        // 获取节点数据：节点索引、节点在 Basic 坐标系下的坐标、节点位移输出坐标系 ID
        ModelInfoQuery::QueryNode (nodeService, element.m_g, 4, _OUT pointSet, _OUT nodeIndexSet, _OUT dispCoordIdSet);

        // 获取单元自由度对应的全局自由度索引
        assembly.GetGlobalDOFId (nodeIndexSet, dofSet, _OUT dofGlobalIndexSet);

        // 计算单元坐标系的变换矩阵和单元坐标系下的节点坐标
        CoordTrans (pointSet, _OUT T, _OUT localCoord);

        // 计算单元的材料坐标系变换矩阵：单元坐标系到材料坐标系
        propertyService->FindById (element.m_pId, _OUT property);
        auto matTypes = materialSummaryService->GetMaterialType (property.m_mId);
        if (std::find (matTypes.begin (), matTypes.end (), MaterialTypes::MATTYPE_3DAnisotropic) != matTypes.end ())
        {
            if (-2 == property.m_mCId)
            {
                ComputeInternalCoordTrans (pointSet, _OUT TS);
            }
            SolidMaterialParser::ComputeMaterialCoordTrans (property.m_mCId, T.Transpose (), TS, _OUT TMat);
        }

        // 计算单元坐标系材料弹性矩阵
        Matrixd D (6, 6);
        SolidMaterialParser::ComputeElasticMatrix (dbServiceFactory, property.m_mId, TMat, _OUT D, _OUT rho);

        Matrixd TEle (12, 12);
        for (std::size_t iNode{0}; iNode < 4; ++iNode)
        {
            const std::size_t i{iNode * 3};
            TEle.BlockFill (i, i, T);
            if (dispCoordIdSet[iNode] > 0)
            {
                coordProcessor.GetCoordTrans (dispCoordIdSet[iNode], pointSet[iNode], _OUT TGSet[iNode]);
                TEle.BlockFill (i, i, TGSet[iNode] * T);
            }
        }

        ComputeStiffness (localCoord, D, _OUT ke);
        ke = TEle * ke * TEle.Transpose ();
        assembly.AccumulateStiffness (dofGlobalIndexSet, ke);

        const auto needCoupledMass = assembly.NeedCoupleMass ();
        if (needCoupledMass)
        {
            mass = ComputeCoupledMass (localCoord, rho);
            if (std::any_of (dispCoordIdSet.begin (), dispCoordIdSet.end (), [] (Id id) {return id > 0;}))
            {
                mass = TEle * mass * TEle.Transpose ();
            }
        }
        else
        {
           ComputeLumpedMass (localCoord, rho, _OUT mass);
        }
        assembly.AccumulateCoupledMass (dofGlobalIndexSet, mass);
    }
}

void Tetra4ElementCalculator::ComputeElementStrainStress
(
const SG::Algebra::Matrixd& localCoord,
const SG::Algebra::Matrixd& D,
const SG::Algebra::Matrixd& elementDisplacement,
const SG::Algebra::Matrixd& TMat,
_OUT SG::Algebra::Matrixd&  elementStrain,
_OUT SG::Algebra::Matrixd&  elementStress
)
{
    // 构建 B 矩阵
    Matrixd B;
    Real V;
    CalcBMatrix (localCoord, _OUT B, _OUT V);

    elementStrain = B * elementDisplacement;
    elementStress = D * elementStrain;

    Matrixd TStrain (6, 6);
    TenserTranform::StrainTrans (TMat, _OUT TStrain);
    elementStrain = TStrain * elementStrain;
    Matrixd TStress (6, 6);
    TenserTranform::StressTrans (TMat, _OUT TStress);
    elementStress = TStress * elementStress;
}

void Tetra4ElementCalculator::DataRecovery
(
const SG::DBManager::DBServiceFactorySP             dbServiceFactory,
int                                                 subcaseId,
SG::DataStructure::Common::Real                     time,
const std::vector<SG::DataStructure::Common::Real>& globalDisplacement
)
{
    const auto& assembly       = AssemblyEigen::GetInstance ();
    const auto& coordProcessor = CoordinateProcessor::GetInstance ();

    static const DofSet dofSet{DOFDirections::DOF_Ux, DOFDirections::DOF_Uy, DOFDirections::DOF_Uz};

    CTETRA4Data              element;
    std::vector<std::size_t> nodeIndexSet (4);
    std::vector<int>         dofGlobalIndexSet (12);
    PointArray               pointSet (4);
    IdArray                  dispCoordIdSet (4);
    PSOLIDData               property;

    Matrixd T (3, 3);
    Matrixd localCoord (3, 4);
    Matrixd TS (3, 3);
    Matrixd TMat (3, 3);
    Real    rho;
    Matrixd eleDisplacement (12, 1);
    Matrixd eleStrain (6, 1);
    Matrixd eleStress (6, 1);

    const auto& nodeService     = dbServiceFactory->get<INodeService> ();
    const auto& eleService      = dbServiceFactory->get<IElementService> ();
    const auto& propertyService = dbServiceFactory->get<IPSOLIDService> ();
    const auto  eleSize         = eleService->GetElementSize (m_type);
    TETRAStrainData eleStrainData;
    std::vector<TETRAStrainData> eleStrainDataSet (eleSize);
    TETRAStressData eleStressData;
    std::vector<TETRAStressData> eleStressDataSet (eleSize);

    for (std::size_t iEle{0}; iEle < eleSize; ++iEle)
    {
        eleService->FindByIndex (iEle, _OUT element);

        // 获取节点数据：节点索引、节点在 Basic 坐标系下的坐标、节点位移输出坐标系 ID
        ModelInfoQuery::QueryNode (nodeService, element.m_g, 4, _OUT pointSet, _OUT nodeIndexSet, _OUT dispCoordIdSet);

        // 获取单元自由度对应的全局自由度索引
        assembly.GetGlobalDOFId (nodeIndexSet, dofSet, _OUT dofGlobalIndexSet);

        // 计算单元坐标系的变换矩阵和单元坐标系下的节点坐标
        CoordTrans (pointSet, _OUT T, _OUT localCoord);

        // 计算单元的材料坐标系变换矩阵：单元坐标系到材料坐标系
        propertyService->FindById (element.m_pId, _OUT property);
        if (-2 == property.m_mCId)
        {
            ComputeInternalCoordTrans (pointSet, _OUT TS);
        }
        SolidMaterialParser::ComputeMaterialCoordTrans (property.m_mCId, T.Transpose (), TS, _OUT TMat);

        // 计算单元坐标系材料弹性矩阵
        Matrixd D (6, 6);
        SolidMaterialParser::ComputeElasticMatrix (dbServiceFactory, property.m_mId, TMat, _OUT D, _OUT rho);

        // 从总体位移向量中提取单元位移
        for (std::size_t iDof{0}; iDof < 12; ++iDof)
        {
            eleDisplacement (iDof, 0) = globalDisplacement[dofGlobalIndexSet[iDof]];
        }

        // 单元位移变换到单元坐标系
        Matrixd TEle (12, 12);
        for (std::size_t iNode{0}; iNode < 4; ++iNode)
        {
            const std::size_t i{iNode * 3};
            TEle.BlockFill (i, i, T);
            if (dispCoordIdSet[iNode] > 0)
            {
                Matrixd TG (3, 3);
                coordProcessor.GetCoordTrans (dispCoordIdSet[iNode], pointSet[iNode], _OUT TG);
                TEle.BlockFill (i, i, TG * T);
            }
        }
        eleDisplacement = TEle.Transpose () * eleDisplacement;

        // 计算单元应力应变
        ComputeElementStrainStress (localCoord, D, eleDisplacement, TMat, _OUT eleStrain, _OUT eleStress);

        // 转换为存储的单元应变数据结构
        eleStrainData.m_eId     = element.m_id;
        eleStrainData.m_CID     = property.m_mCId;
        eleStrainData.SetCType (property.m_stress);
        eleStrainData.m_NODEF   = 4;
        eleStrainData.m_GRID[0] = 0;
        eleStrainData.m_X[0]    = eleStrain (0, 0);
        eleStrainData.m_Y[0]    = eleStrain (1, 0);
        eleStrainData.m_Z[0]    = eleStrain (2, 0);
        eleStrainData.m_TXY[0]  = eleStrain (3, 0);
        eleStrainData.m_TYZ[0]  = eleStrain (4, 0);
        eleStrainData.m_TZX[0]  = eleStrain (5, 0);
        for (std::size_t iNode{1}; iNode < 5; ++iNode)
        {
            eleStrainData.m_GRID[iNode] = element.m_g[iNode-1];
            eleStrainData.m_X[iNode]    = eleStrain (0, 0);
            eleStrainData.m_Y[iNode]    = eleStrain (1, 0);
            eleStrainData.m_Z[iNode]    = eleStrain (2, 0);
            eleStrainData.m_TXY[iNode]  = eleStrain (3, 0);
            eleStrainData.m_TYZ[iNode]  = eleStrain (4, 0);
            eleStrainData.m_TZX[iNode]  = eleStrain (5, 0);
        }
        eleStrainDataSet[iEle] = eleStrainData;

        // 转换为存储的单元应力数据结构
        eleStressData.m_eId     = element.m_id;
        eleStressData.m_CID     = property.m_mCId;
        eleStressData.SetCType (property.m_stress);
        eleStressData.m_NODEF   = 4;
        eleStressData.m_GRID[0] = 0;
        eleStressData.m_X[0]    = eleStress (0, 0);
        eleStressData.m_Y[0]    = eleStress (1, 0);
        eleStressData.m_Z[0]    = eleStress (2, 0);
        eleStressData.m_TXY[0]  = eleStress (3, 0);
        eleStressData.m_TYZ[0]  = eleStress (4, 0);
        eleStressData.m_TZX[0]  = eleStress (5, 0);
        for (std::size_t iNode{1}; iNode < 5; ++iNode)
        {
            eleStressData.m_GRID[iNode] = element.m_g[iNode-1];
            eleStressData.m_X[iNode]    = eleStress (0, 0);
            eleStressData.m_Y[iNode]    = eleStress (1, 0);
            eleStressData.m_Z[iNode]    = eleStress (2, 0);
            eleStressData.m_TXY[iNode]  = eleStress (3, 0);
            eleStressData.m_TYZ[iNode]  = eleStress (4, 0);
            eleStressData.m_TZX[iNode]  = eleStress (5, 0);
        }
        eleStressDataSet[iEle] = eleStressData;
    }

    dbServiceFactory->get<ITETRAStrainService> ()->AddBySubcaseId (static_cast<Id> (subcaseId), eleStrainDataSet);
    dbServiceFactory->get<ITETRAStressService> ()->AddBySubcaseId (static_cast<Id> (subcaseId), eleStressDataSet);
}

void Tetra4ElementCalculator::ComputeStressStiffness (const DBServiceFactorySP dbServiceFactory, int staticSubId)
{
    auto&       assembly       = AssemblyEigen::GetInstance ();
    const auto& coordProcessor = CoordinateProcessor::GetInstance ();

    static const DofSet dofSet{DOFDirections::DOF_Ux, DOFDirections::DOF_Uy, DOFDirections::DOF_Uz};

    CTETRA4Data              element;
    std::vector<std::size_t> nodeIndexSet (4);
    std::vector<int>         dofGlobalIndexSet (12);
    PointArray               pointSet (4);
    IdArray                  dispCoordIdSet (4);
    PSOLIDData               property;

    Matrixd T (3, 3);
    Matrixd localCoord (3, 4);
    Matrixd TS (3, 3);
    Matrixd kSigma (12, 12);
    Matrixd eleMatCoordStress (6, 1);
    Matrixd eleEleCoordStress (6, 1);

    const auto& nodeService     = dbServiceFactory->get<INodeService> ();
    const auto& eleService      = dbServiceFactory->get<IElementService> ();
    const auto& propertyService = dbServiceFactory->get<IPSOLIDService> ();
    const auto  eleSize         = eleService->GetElementSize (m_type);

    //从数据库提取积分点应力结果
    std::vector<TETRAStressData> eleStressDataSet;
    dbServiceFactory->get<ITETRAStressService> ()->Find ({SOLType::SOL_SESTATIC, static_cast<Id> (staticSubId)}, _OUT eleStressDataSet);
    
    for (std::size_t iEle{0}; iEle < eleSize; ++iEle)
    {
        eleService->FindByIndex (iEle, _OUT element);

        // 获取节点数据：节点索引、节点在 Basic 坐标系下的坐标、节点位移输出坐标系 ID
        ModelInfoQuery::QueryNode (nodeService, element.m_g, 4, pointSet, nodeIndexSet, dispCoordIdSet);

        // 获取单元自由度对应的全局自由度索引
        assembly.GetGlobalDOFId (nodeIndexSet, dofSet, _OUT dofGlobalIndexSet);

        // 计算单元坐标系的变换矩阵和单元坐标系下的节点坐标
        CoordTrans (pointSet, _OUT T, _OUT localCoord);

        const auto& eleStressData = eleStressDataSet[iEle];
        //单元材料坐标系下积分点应力
        eleMatCoordStress (0, 0) = eleStressData.m_X[0];
        eleMatCoordStress (1, 0) = eleStressData.m_Y[0];
        eleMatCoordStress (2, 0) = eleStressData.m_Z[0];
        eleMatCoordStress (3, 0) = eleStressData.m_TXY[0];
        eleMatCoordStress (4, 0) = eleStressData.m_TYZ[0];
        eleMatCoordStress (5, 0) = eleStressData.m_TZX[0];

        //应力参考坐标系由材料坐标系转换至单元坐标系
        //单元坐标系下应力是DataRecovery计算材料坐标系下应力的中间量，此处重复计算后续可优化

        // 材料弹性矩阵
        propertyService->FindById (element.m_pId, _OUT property);
        //此处计算流程与单元应力转换至材料局部坐标应力相同
        if (-2 == property.m_mCId)
        {
            ComputeInternalCoordTrans (pointSet, _OUT TS);
        }
        Matrixd TMat (3, 3);
        SolidMaterialParser::ComputeMaterialCoordTrans (property.m_mCId, T.Transpose (), TS, _OUT TMat);

        Matrixd TStress (6, 6);
        TenserTranform::StrainTrans (TMat, _OUT TStress);
        //涉及求逆，计算量较大且此处为重复计算，可优化
        eleEleCoordStress = TStress.Transpose () * eleMatCoordStress;

        ComputeStressStiffness (localCoord, eleEleCoordStress, _OUT kSigma);

        // 单元刚度矩阵变换到位移输出坐标系
        Matrixd TEle (12, 12);
        for (std::size_t iNode{0}; iNode < 4; ++iNode)
        {
            const std::size_t i{iNode * 3};
            TEle.BlockFill (i, i, T);
            if (dispCoordIdSet[iNode] > 0)
            {
                Matrixd TG (3, 3);
                coordProcessor.GetCoordTrans (dispCoordIdSet[iNode], pointSet[iNode], _OUT TG);
                TEle.BlockFill (i, i, TG * T);
            }
        }
        kSigma = TEle * kSigma * TEle.Transpose ();

        assembly.AccumulateStressStiffness (dofGlobalIndexSet, kSigma);
    }
}